home *** CD-ROM | disk | FTP | other *** search
/ Aminet 2 / Aminet AMIGA CDROM (1994)(Walnut Creek)[Feb 1994][W.O. 44790-1].iso / Aminet / util / misc / Fudgit233.lha / Source / src / fits.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-12-16  |  13.2 KB  |  430 lines
  1. #include <math.h>
  2. #include <stdio.h>
  3. #include <string.h>
  4. #ifndef NOUNISTD_H
  5. #include <unistd.h>
  6. #endif
  7.  
  8. #include "symbol.h"
  9. #include "code.h"
  10. #include "math.tab.h"
  11. #include "fudgit.h"
  12. #include "setshow.h"
  13. #include "head.h"
  14.  
  15. static double *Pvec[4];
  16. static void yfit(double **der, double *par, double *vec, int ndata);
  17. static void fexpo(double *x, double *a, double *y, double **dyda, int ndata);
  18. static void fcosine(double *X, double **du, int ndata);
  19. static void fsine(double *X, double **du, int ndata);
  20. static void fgauss(double *x, double *a, double *y, double **dyda, int ndata);
  21. static void fpoly(double *X, double **du, int ndata);
  22. static void fleg(double *X, double **du, int ndata);
  23.  
  24. extern double *Ft_X2;
  25. extern double *Ft_Param;
  26. extern void Ft_mathuser(void);
  27.  
  28.  
  29. extern void Ft_fit (double *x, double *y, int ndata, double *sig, int mwt, double *a, double *b, double *siga, double *sigb, double *chi2, double *q);
  30. extern void Ft_pearsn (double *x, double *y, int n, double *r, double *prob, double *z);
  31. extern void Ft_medfit (double *x, double *y, int ndata, double *a, double *b, double *abdev);
  32. extern int Ft_svdvar (double **v, int ma, double *w, double **cvm);
  33.  
  34. int Ft_fits(int argc, char **argv, int ndata)
  35. {
  36.     int i, j;
  37.     Symbol *sym[4];
  38.     char name[TOKENSIZE+10];
  39.  
  40.     extern int Ft_svdfit(double **coef, double *y, double *sig, int ndata, double *a, int ma, double **v, double *w, double *chisq);
  41.     extern int Ft_lfit(double **coef, double *y, double *sig, int ndata, double *a, int ma, int *lista, int mfit, double **covar, double *chisq);
  42.     extern int Ft_mrqmin(double **coef, double *f, double *y, double *sig, int ndata, double *a, int ma, int *lista, int mfit, double **covar, double **alpha, double *chisq, double *alamda);
  43.     extern Symbol *Ft_lookup(char *s);
  44.    
  45. /* printf("CHECK FOR STUPID MISTAKES\n"); */
  46.     if (ndata == 0) {
  47.         fputs("No data to fit!\n", stderr);
  48.         return(ERRR);
  49.     }
  50.     else if ((int) *Ft_Param == 0) {
  51.         fputs("You must set parameters first!\n", stderr);
  52.         return(ERRR);
  53.     }
  54.     else if (Ft_Methi == NONE) {
  55.         fputs("You must set method first!\n", stderr);
  56.         return(ERRR);
  57.     }
  58.     else if (Ft_Funci == NONE) {
  59.         fputs("You must set function first!\n", stderr);
  60.         return(ERRR);
  61.     }
  62.     else if ((Ft_Methi == LS_FIT || Ft_Methi == ML_FIT) && Ft_Mlist == 0) {
  63.         fprintf(stderr, "You must 'adjust' before fitting with '%s'.\n",
  64.         Ft_Method[Ft_Methi].name);
  65.         return(ERRR);
  66.     }
  67. /* printf("LOOKUP-INSTALL FIT VARIABLES\n"); */
  68.     for (i=0;i<argc-1;i++) {
  69.         if ((sym[i] = Ft_lookup(argv[i+1])) == 0) {
  70.             fprintf(stderr, "%s: No such variable.\n", argv[i+1]);
  71.             return(ERRR);
  72.         }
  73.         else if (sym[i]->type != VEC) {
  74.             fprintf(stderr, "%s: Not a vector.\n", argv[i+1]);
  75.             return(ERRR);
  76.         }
  77.         Pvec[i] = sym[i]->u.vec;
  78.     }
  79. /* printf("make fit vector\n"); */
  80.     sprintf(name, "%sFIT", sym[1]->name);
  81.     if ((sym[3] = Ft_lookup(name)) == 0) {
  82.         sym[3] = Ft_install(name, VEC, Ft_Samples);
  83.     }
  84.     Pvec[3] = sym[3]->u.vec;
  85.     for (i=1;i<= (int) *Ft_Param;i++) {
  86.         sprintf(name, "D%sFITD%d", sym[1]->name, i);
  87.         if ((sym[3] = Ft_lookup(name)) == 0) {
  88.             sym[3] = Ft_install(name, VEC, Ft_Samples);
  89.         }
  90.         Ft_Mparxsamp[i] = sym[3]->u.vec;
  91.     }
  92.    
  93. /* printf("PREPARE PARTIAL DERIVATIVE MATRIX\n"); */
  94.     /* here we start --- going Ft_Functionwise   */
  95.     if (Ft_Funci == POLY) {
  96.         fpoly(Pvec[0], Ft_Mparxsamp, ndata);
  97.     }
  98.     else if (Ft_Funci == SINE) {
  99.         fsine(Pvec[0], Ft_Mparxsamp, ndata);
  100.     }
  101.     else if (Ft_Funci == COSINE) {
  102.         fcosine(Pvec[0], Ft_Mparxsamp, ndata);
  103.     }
  104.     else if (Ft_Funci == LEG) {
  105.         fleg(Pvec[0], Ft_Mparxsamp, ndata);
  106.     }
  107.     else if (Ft_Funci == GAUSS) {
  108.         if ((int) *Ft_Param%3) {
  109.             fputs("Number of parameters not a multiple of 3.\n", stderr);
  110.             return(ERRR);
  111.         }
  112.         if (Ft_Methi != ML_FIT) {
  113.             fprintf(stderr, "'%s' can only be used with '%s'.\n",
  114.             Ft_Function[Ft_Funci].name, Ft_Method[ML_FIT].name);
  115.             return(ERRR);
  116.         }
  117.         fgauss(Pvec[0], Ft_A, Pvec[3], Ft_Mparxsamp, ndata);
  118.     }
  119.     else if (Ft_Funci == EXPO) {
  120.         if ((int) *Ft_Param%2) {
  121.             fputs("Number of parameters not a multiple of 2.\n", stderr);
  122.             return(ERRR);
  123.         }
  124.         if (Ft_Methi != ML_FIT) {
  125.             fprintf(stderr, "'%s' can only be used with '%s'.\n",
  126.             Ft_Function[Ft_Funci].name, Ft_Method[ML_FIT].name);
  127.             return(ERRR);
  128.         }
  129.         fexpo(Pvec[0], Ft_A, Pvec[3], Ft_Mparxsamp, ndata);
  130.     }
  131.     else if (Ft_Funci == STRL) {
  132. /* printf("IF STRAIGHT LINE DO IT NOW\n"); */
  133.         if (Ft_Methi == SVD_FIT || Ft_Methi == LS_FIT) {
  134.             fprintf(stderr, "Use a polynomial for '%s'!\n",
  135.             Ft_Method[Ft_Methi].name);
  136.             return(ERRR);
  137.         }
  138.         else if (Ft_Methi == LS_REG) {
  139.             int dy = 1;
  140.  
  141.             if (argc == 3) {
  142.                 dy = 0;
  143.                 Pvec[2] = 0;
  144.             }
  145.             if ((int) *Ft_Param != 2) {
  146.                 fprintf(stderr, "2 parameters are needed for '%s'.\n",
  147.                 Ft_Method[LS_REG].name);
  148.                 return(ERRR);
  149.             }
  150.             fpoly(Pvec[0], Ft_Mparxsamp, ndata);
  151.             Ft_fit(Pvec[0], Pvec[1], ndata, Pvec[2], dy,
  152.             Ft_A+1, Ft_A+2, Ft_DA+1, Ft_DA+2, Ft_X2, &Ft_Q);
  153.             yfit(Ft_Mparxsamp, Ft_A, Pvec[3], ndata);
  154.             Ft_pearsn(Pvec[0], Pvec[1], ndata, Ft_Cortest, Ft_Cortest+1,
  155.             Ft_Cortest+2);
  156.         }
  157.         else if (Ft_Methi == LA_REG) {
  158.             if ((int) *Ft_Param != 2) {
  159.                 fprintf(stderr, "2 parameters are needed for '%s'.\n",
  160.                 Ft_Method[LA_REG].name);
  161.                 return(ERRR);
  162.             }
  163.             if (argc == 4) {
  164.                 fprintf(stderr, "Warning: Errors are not used in '%s'.\n",
  165.                 Ft_Method[LA_REG].name);
  166.             }
  167.             fpoly(Pvec[0], Ft_Mparxsamp, ndata);
  168.             Ft_medfit(Pvec[0], Pvec[1], ndata, Ft_A+1, Ft_A+2, Ft_X2);
  169.             yfit(Ft_Mparxsamp, Ft_A, Pvec[3], ndata);
  170.             Ft_pearsn(Pvec[0], Pvec[1], ndata, Ft_Cortest,
  171.             Ft_Cortest+1, Ft_Cortest+2);
  172.         }
  173.         else {
  174.             fprintf(stderr, "Ft_Function '%s' is not supported under '%s'.\n",
  175.             Ft_Function[Ft_Funci].name, Ft_Method[Ft_Methi].name);
  176.             return(ERRR);
  177.         }
  178.         return(0);
  179.     }
  180.     else if (Ft_Funci == USER && Ft_Methi != ML_FIT) {
  181.         Ft_mathuser();
  182.     }
  183.  
  184. /* printf("WHICH OTHER METHOD\n"); */
  185.     /* Now going method wise  */
  186.     if (Ft_Methi == SVD_FIT) {
  187.         if (Ft_svdfit(Ft_Mparxsamp, Pvec[1], Pvec[2],
  188.         ndata, Ft_A, (int) *Ft_Param, Ft_M2parxpar, Ft_Mfparx1, Ft_X2) != ERRR) {
  189.             Ft_svdvar(Ft_M2parxpar, (int) *Ft_Param, Ft_Mfparx1, Ft_M1parxpar);
  190.             yfit(Ft_Mparxsamp, Ft_A, Pvec[3], ndata);
  191.             for (i=1;i<= (int) *Ft_Param;i++) {
  192.                 Ft_DA[i] = sqrt(Ft_M1parxpar[i][i]);
  193.             }
  194.         }
  195.     }
  196.     else if (Ft_Methi == LS_FIT) {
  197.         if (Ft_lfit(Ft_Mparxsamp, Pvec[1], Pvec[2], ndata,
  198.         Ft_A, (int) *Ft_Param, Ft_Miparx1, Ft_Mlist, Ft_M1parxpar, Ft_X2) != ERRR) {
  199.             yfit(Ft_Mparxsamp, Ft_A, Pvec[3], ndata);
  200.             for (i=1;i<= (int) *Ft_Param;i++) {
  201.                 Ft_DA[i] = sqrt(Ft_M1parxpar[i][i]);
  202.             }
  203.         }
  204.     }
  205.     else if (Ft_Methi == ML_FIT) {
  206.         int same = 0;
  207.         double ox2 = 0.0;
  208.         int iter = (Ft_Iter > 0? Ft_Iter : -Ft_Iter);
  209.  
  210.         Ft_Q = -1;
  211.         if (Ft_mrqmin(Ft_Mparxsamp, Pvec[3], Pvec[1], Pvec[2], ndata, Ft_A,
  212.         (int) *Ft_Param, Ft_Miparx1, Ft_Mlist, Ft_M1parxpar,
  213.         Ft_M2parxpar, Ft_X2, &Ft_Q) == ERRR) {
  214.             return(ERRR);
  215.         }
  216.         for (i=1; i<= iter;i++) {
  217.             if (Ft_mrqmin(Ft_Mparxsamp, Pvec[3], Pvec[1], Pvec[2], ndata, Ft_A,
  218.             (int) *Ft_Param, Ft_Miparx1, Ft_Mlist, Ft_M1parxpar, Ft_M2parxpar,
  219.             Ft_X2, &Ft_Q) == ERRR) {
  220.                 break;
  221.             }
  222.             fprintf(stderr, "~~~~~~~~~~ Iteration: %02d ~~~~~~~~~~~~\t", i);
  223.             fprintf(stderr, "%s: ", "Chi^2");
  224.             fprintf(stderr, Ft_Format, *Ft_X2);
  225.             fprintf(stderr, "\n");
  226.             ox2 = *Ft_X2 -ox2;
  227.             fprintf(stderr, "%45s: ", "DChi^2");
  228.             fprintf(stderr, Ft_Format, ox2);
  229.             fputc('\n', stderr);
  230.             for (j=1;j <= (int) *Ft_Param;j++) {
  231.                 char buffer[128];
  232.                 sprintf(buffer, "%s[%d]", Ft_Pname, j);
  233.                 fprintf(stderr, "%20s: ", buffer);
  234.                 fprintf(stderr, Ft_Format, Ft_A[j]);
  235.                 fputc('\n', stderr);
  236.             }
  237.             if (Ft_Iter > 0) { /* check for convergence if Iter > 0 */
  238.                 if (ox2 == (double)0.0)
  239.                     same++;
  240.                 else
  241.                     same = 0;
  242.                 if (same == 2) {
  243.                     fprintf(stderr, "Convergence after %d iterations.\n", i);
  244.                     break;
  245.                 }
  246.                 ox2 = *Ft_X2;
  247.             }
  248.         }
  249.         Ft_Q = 0.0;
  250.         if (Ft_mrqmin(Ft_Mparxsamp, Pvec[3], Pvec[1], Pvec[2], ndata, Ft_A,
  251.         (int) *Ft_Param, Ft_Miparx1, Ft_Mlist, Ft_M1parxpar,
  252.         Ft_M2parxpar, Ft_X2, &Ft_Q) == ERRR) {
  253.             return(ERRR);
  254.         }
  255.         for (i=1;i<= (int) *Ft_Param;i++) {
  256.             Ft_DA[i] = sqrt(Ft_M1parxpar[i][i]);
  257.         }
  258.     }
  259.     else {
  260.         fprintf(stderr, "Ft_Function '%s' is not supported under '%s'.\n",
  261.         Ft_Function[Ft_Funci].name, Ft_Method[Ft_Methi].name);
  262.         return(ERRR);
  263.     }
  264.     return(0);
  265. }
  266.  
  267. void Ft_fml_fit(int ndata)
  268. {
  269.     extern void Ft_mathuser(void);
  270.  
  271.     switch(Ft_Funci) {
  272.         case USER:
  273.             Ft_mathuser();
  274.             return;
  275.         case GAUSS:
  276.             fgauss(Pvec[0], Ft_A, Pvec[3], Ft_Mparxsamp, ndata);
  277.             return;
  278.         case EXPO:
  279.             fexpo(Pvec[0], Ft_A, Pvec[3], Ft_Mparxsamp, ndata);
  280.             return;
  281.     }
  282. }
  283.  
  284. /* linear */
  285. static void fleg(double *X, double **du, int ndata)
  286. {
  287.     int i, j;
  288.     double twox,f2,f1,d;
  289.  
  290.     if (Ft_Methi != LS_FIT && Ft_Methi != SVD_FIT) {
  291.         fprintf(stderr,
  292.         "Legendre series not supported under '%s' method.\n",
  293.         Ft_Method[Ft_Methi].name);
  294.         Ft_catcher(ERRR);
  295.     }
  296.     for (i=1; i<= ndata; i++) {
  297.         du[1][i]=1.0;
  298.         du[2][i]=X[i];
  299.         if ((int) *Ft_Param > 2) {
  300.             twox=2.0*X[i];
  301.             f2=X[i];
  302.             d=1.0;
  303.             for (j=3;j<=(int) *Ft_Param;j++) {
  304.                 f1=d;
  305.                 d += 1.0;
  306.                 f2 += twox;
  307.                 du[j][i]=(f2*du[j-1][i]-f1*du[j-2][i])/d;
  308.             }
  309.         }
  310.     }
  311.     return;
  312. }
  313.  
  314. /* linear */
  315. static void fpoly(double *X, double **du, int ndata)
  316. {
  317.     int i, j;
  318.  
  319.     if (Ft_Methi == ML_FIT) {
  320.         fprintf(stderr,
  321.         "Power series not supported under '%s' method.\n",
  322.         Ft_Method[Ft_Methi].name);
  323.         Ft_catcher(ERRR);
  324.     }
  325.     for (i=1;i<=ndata;i++) {
  326.         du[1][i]=1.0;
  327.         for (j=2;j<=(int) *Ft_Param;j++) du[j][i]=du[j-1][i]*X[i];
  328.     }
  329.     return;
  330. }
  331.  
  332. static void fgauss(double *x, double *a, double *y, double **dyda, int ndata)
  333. /* non-linear  */
  334. {
  335.     int i, j;
  336.     double fac, ex, arg;
  337.  
  338.     for (j=1;j<= ndata;j++) {
  339.         y[j]=0.0;
  340.         for (i=1;i<=(int) *Ft_Param-1;i+=3) {
  341.             arg=(x[j]-a[i+1])/a[i+2];
  342.             ex=exp(-arg*arg);
  343.             fac=a[i]*ex*2.0*arg;
  344.             y[j] += a[i]*ex;
  345.             dyda[i][j]=ex;
  346.             dyda[i+1][j]=fac/a[i+2];
  347.             dyda[i+2][j]=fac*arg/a[i+2];
  348.         }
  349.     }
  350.     return;
  351. }
  352.  
  353. static void fsine(double *X, double **du, int ndata)  /* linear */
  354.                 
  355.           
  356. {
  357.     int i, j;
  358.  
  359.     if (Ft_Methi != LS_FIT  && Ft_Methi != SVD_FIT) {
  360.         fprintf(stderr,
  361.         "Sine series not supported under '%s' method.\n",
  362.         Ft_Method[Ft_Methi].name);
  363.         Ft_catcher(ERRR);
  364.     }
  365.     for (i=1;i<= ndata;i++) {
  366.         for (j=1;j<= (int) *Ft_Param;j++) {
  367.             du[j][i] = sin(j*X[i]);
  368.         }
  369.     }
  370.     return;
  371. }
  372.  
  373. static void fcosine(double *X, double **du, int ndata)  /* linear */
  374.                 
  375.           
  376. {
  377.     int i, j;
  378.  
  379.     if (Ft_Methi != LS_FIT  && Ft_Methi != SVD_FIT) {
  380.         fprintf(stderr,
  381.         "Sine series not supported under '%s' method.\n",
  382.         Ft_Method[Ft_Methi].name);
  383.         Ft_catcher(ERRR);
  384.     }
  385.     for (i=1;i<= ndata;i++) {
  386.         du[1][i] = 1.0;
  387.         for (j=2; j<= (int) *Ft_Param;j++) {
  388.             du[j][i] = cos((j-1)*X[i]);
  389.         }
  390.     }
  391.     return;
  392. }
  393.  
  394. static void fexpo(double *x, double *a, double *y, double **dyda, int ndata)   /* non-linear  */
  395.                           
  396.           
  397. {
  398.     int i, j;
  399.     double ex, arg;
  400.  
  401.     for (j=1;j <= ndata;j++) {
  402.         y[j] = 0.0;
  403.         for (i=1;i<=(int) *Ft_Param-1;i+=2) {
  404.             arg = (x[j] * a[i+1]);
  405.             ex = exp(arg);
  406.             y[j] += a[i] * ex;
  407.             dyda[i][j] = ex;
  408.             dyda[i+1][j] = a[i] * a[i+1] * ex;
  409.         }
  410.     }
  411.     return;
  412. }
  413.  
  414. /* rebuild the function from the partial(total) linear derivatives  */
  415. /* not valid for non-linear fits  */
  416. static void yfit(double **der, double *par, double *vec, int ndata)
  417. {
  418.     int i, j;
  419.  
  420.     for (j=1;j<=ndata;j++) {
  421.         vec[j] = 0.0;
  422.     }
  423.     for (i=1;i<=(int) *Ft_Param;i++) {
  424.         for (j=1;j<=ndata;j++) {
  425.             vec[j] += par[i]*der[i][j];
  426.         }
  427.     }
  428.     return;
  429. }
  430.